Fijamos la semilla para poder reproducir los experimentos. También se fija el numero de CPU’s a utilizar.
set.seed(42)
options(mc.cores = 24)
Se importan las librerÃas a utilizar a lo largo de la notebook:
# install.packages(pacman)
# install.packages("https://cran.r-project.org/src/contrib/rstan_2.21.2.tar.gz",repos = NULL,type="source")
# sudo apt-get install libglpk-dev
library(pacman)
p_load(tidyverse, tidymodels, rsample, rstan, shinystan, rstanarm, devtools)
source('../src/dataset.R')
source('../src/plot.R')
source('../src/model.R')
Palmer Penguins
dataset <- load_dataset() %>% mutate_if(is.character, as.factor)
Rows: 344 Columns: 8
── Column specification ───────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (3): species, island, sex
dbl (5): bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g, year
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dataset %>% glimpse()
Rows: 344
Columns: 8
$ species <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adelie,…
$ island <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torgersen, Torgersen, Torgersen, To…
$ bill_length_mm <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.1, 42.0, 37.8, 37.8, 41.1, 38.…
$ bill_depth_mm <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.1, 20.2, 17.1, 17.3, 17.6, 21.…
$ flipper_length_mm <dbl> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, 186, 180, 182, 191, 198, 185, …
$ body_mass_g <dbl> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 3475, 4250, 3300, 3700, 3200, 380…
$ sex <fct> male, female, female, NA, female, male, female, male, NA, NA, NA, NA, female, m…
$ year <dbl> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2…
Se enumeran y describen breve-mente cada variable que forma parte del dataset:
Variables Numéricas:
Variables Categóricas:
A continuación veamos los posibles valores de las variables categóricas:
show_values(dataset %>% select_if(negate(is.numeric)))
_____________
species n
=============
Adelie 152
Chinstrap 68
Gentoo 124
¯¯¯¯¯¯¯¯¯¯¯¯¯
_____________
island n
=============
Biscoe 168
Dream 124
Torgersen 52
¯¯¯¯¯¯¯¯¯¯¯¯¯
__________
sex n
==========
female 165
male 168
NA 11
¯¯¯¯¯¯¯¯¯¯
missings_summary(dataset)
hist_plots(dataset)
Observaciones
box_plots(dataset)
Observaciones
No se registran valores mas extremos que el mÃnimo y máximo valor en cada variables. Es decir que no encontramos outliers.
outliers(dataset, column='bill_length_mm')
$inf
[1] 32.1
$sup
[1] 59.6
outliers(dataset, column='bill_length_mm')
$inf
[1] 32.1
$sup
[1] 59.6
outliers(dataset, column='bill_depth_mm')
$inf
[1] 13.1
$sup
[1] 21.5
outliers(dataset, column='flipper_length_mm')
$inf
[1] 172
$sup
[1] 231
outliers(dataset, column='body_mass_g')
$inf
[1] 2700
$sup
[1] 6300
outliers(dataset, column='year')
$inf
[1] 2007
$sup
[1] 2009
bar_plots(dataset)
Observaciones
dataset <- dataset %>% drop_na()
missings_summary(dataset)
Warning: attributes are not identical across measure variables;
they will be dropped
corr_plot(dataset %>% dplyr::select(-year))
segmented_pairs_plot(dataset, segment_column='species')
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
train_test <- train_test_split(dataset, train_size = 0.7, shuffle = TRUE)
[1] "Train set size: 233"
[1] "Test set size: 100"
train_set <- train_test[[1]]
test_set <- train_test[[2]]
lineal_model_1 <- lm(
body_mass_g ~ bill_length_mm + bill_depth_mm + flipper_length_mm,
data = train_set
)
bayesion_model_1 <- stan(
model_code = "
data {
int<lower=1> obs_count;
vector<lower=1>[obs_count] x1;
vector<lower=1>[obs_count] x2;
vector<lower=1>[obs_count] x3;
vector[obs_count] y;
}
parameters {
real beta0;
real beta1;
real beta2;
real beta3;
real<lower=0> sigma;
}
model {
// Distribuciones a priori
beta0 ~ normal(0, 8000 /* Fijado por el investigador */ ); // Intercept
beta1 ~ normal(0, 100 /* Fijado por el investigador */);
beta2 ~ normal(0, 100 /* Fijado por el investigador */);
beta3 ~ normal(0, 100 /* Fijado por el investigador */);
sigma ~ exponential(0.1 /* Fijado por el investigador */);
// Likelihood
y ~ normal(beta0 + beta1 * x1 + beta2 * x2 + beta3 * x3, sigma);
}
",
data = list(
obs_count = nrow(train_set),
y = colvalues(train_set, 'body_mass_g'),
x1 = colvalues(train_set, 'bill_length_mm'),
x2 = colvalues(train_set, 'bill_depth_mm'),
x3 = colvalues(train_set, 'flipper_length_mm')
),
chains = 3,
iter = 300,
warmup = 180,
thin = 1
)
starting worker pid=2744 on localhost:11086 at 19:38:44.884
starting worker pid=2781 on localhost:11086 at 19:38:44.993
starting worker pid=2818 on localhost:11086 at 19:38:45.103
SAMPLING FOR MODEL 'ff7fc5db0fd7100c79e067b72d83d19f' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 3.7e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.37 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 300 [ 0%] (Warmup)
Chain 1: Iteration: 30 / 300 [ 10%] (Warmup)
Chain 1: Iteration: 60 / 300 [ 20%] (Warmup)
Chain 1: Iteration: 90 / 300 [ 30%] (Warmup)
Chain 1: Iteration: 120 / 300 [ 40%] (Warmup)
Chain 1: Iteration: 150 / 300 [ 50%] (Warmup)
Chain 1: Iteration: 180 / 300 [ 60%] (Warmup)
Chain 1: Iteration: 181 / 300 [ 60%] (Sampling)
Chain 1: Iteration: 210 / 300 [ 70%] (Sampling)
Chain 1: Iteration: 240 / 300 [ 80%] (Sampling)
SAMPLING FOR MODEL 'ff7fc5db0fd7100c79e067b72d83d19f' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 4.7e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.47 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 300 [ 0%] (Warmup)
Chain 2: Iteration: 30 / 300 [ 10%] (Warmup)
Chain 1: Iteration: 270 / 300 [ 90%] (Sampling)
Chain 1: Iteration: 300 / 300 [100%] (Sampling)
Chain 2: Iteration: 60 / 300 [ 20%] (Warmup)
Chain 1:
Chain 1: Elapsed Time: 0.84136 seconds (Warm-up)
Chain 1: 0.717478 seconds (Sampling)
Chain 1: 1.55884 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'ff7fc5db0fd7100c79e067b72d83d19f' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 4.4e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.44 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 300 [ 0%] (Warmup)
Chain 3: Iteration: 30 / 300 [ 10%] (Warmup)
Chain 3: Iteration: 60 / 300 [ 20%] (Warmup)
Chain 2: Iteration: 90 / 300 [ 30%] (Warmup)
Chain 3: Iteration: 90 / 300 [ 30%] (Warmup)
Chain 2: Iteration: 120 / 300 [ 40%] (Warmup)
Chain 3: Iteration: 120 / 300 [ 40%] (Warmup)
Chain 2: Iteration: 150 / 300 [ 50%] (Warmup)
Chain 3: Iteration: 150 / 300 [ 50%] (Warmup)
Chain 2: Iteration: 180 / 300 [ 60%] (Warmup)
Chain 2: Iteration: 181 / 300 [ 60%] (Sampling)
Chain 3: Iteration: 180 / 300 [ 60%] (Warmup)
Chain 3: Iteration: 181 / 300 [ 60%] (Sampling)
Chain 2: Iteration: 210 / 300 [ 70%] (Sampling)
Chain 3: Iteration: 210 / 300 [ 70%] (Sampling)
Chain 2: Iteration: 240 / 300 [ 80%] (Sampling)
Chain 3: Iteration: 240 / 300 [ 80%] (Sampling)
Chain 2: Iteration: 270 / 300 [ 90%] (Sampling)
Chain 3: Iteration: 270 / 300 [ 90%] (Sampling)
Chain 2: Iteration: 300 / 300 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 1.16549 seconds (Warm-up)
Chain 2: 0.466704 seconds (Sampling)
Chain 2: 1.63219 seconds (Total)
Chain 2:
Chain 3: Iteration: 300 / 300 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.802516 seconds (Warm-up)
Chain 3: 0.388908 seconds (Sampling)
Chain 3: 1.19142 seconds (Total)
Chain 3:
params_1 <- c('beta0', 'beta1', 'beta2', 'beta3', 'sigma')
traceplot(bayesion_model_1, pars = params_1, inc_warmup = TRUE)
lm_vs_br_coeficients(lineal_model_1, bayesion_model_1, params_1)
vars_1 <- c('bill_length_mm', 'bill_depth_mm', 'flipper_length_mm')
lm_vs_br_models_validation(lineal_model_1, bayesion_model_1, params_1, vars_1, test_set)
bayesion_predictor_1 <- BayesianRegressionPredictor.from(bayesion_model_1, params_1, vars_1)
plot_compare_fit(
lineal_model_1,
bayesion_predictor_1,
train_set,
label_1='Regresion Lineal',
label_2='Regresion Bayesiana'
)
lineal_model_2 <- lm(
body_mass_g
~ bill_length_mm
+ bill_depth_mm
+ flipper_length_mm
+ sex,
data = train_set
)
Antes que nada transformamos la columna categórica a one-hot encoding:
cat_cols <- c('sex')
train_set_cat <- train_set %>% dummify(cat_cols)
test_set_cat <- test_set %>% dummify(cat_cols)
train_set_cat
Construimos una matriz con todas las variables(X) mas el intercept:
to_model_input <- function(df) {
model.matrix(
body_mass_g
~ bill_length_mm
+ bill_depth_mm
+ flipper_length_mm
+ sex_female
+ sex_male,
data = df
)
}
train_X <- train_set_cat %>% to_model_input()
test_X <- test_set_cat %>% to_model_input()
Definimos y corremos el modelo bayesiano:
bayesion_model_2 <- stan(
model_code = "
data {
int<lower=1> obs_count;
int<lower=1> coef_count;
matrix[obs_count,coef_count] X;
vector[obs_count] y;
}
parameters {
vector[coef_count] beta;
real<lower=0> sigma;
}
model {
// Distribuciones a priori
beta[1] ~ normal(0, 2000); // Intecept
// Variables numericas
beta[2] ~ normal(0, 30);
beta[3] ~ normal(0, 100);
beta[4] ~ normal(0, 100);
// Sexo
beta[5] ~ normal(0, 100);
beta[6] ~ normal(0, 100);
sigma ~ exponential(0.1);
// Likelihood
y ~ normal(X * beta, sigma);
}
",
data = list(
obs_count = dim(train_X)[1],
coef_count = dim(train_X)[2],
y = colvalues(train_set_cat, 'body_mass_g'),
X = train_X
),
chains = 3,
iter = 300,
warmup = 180,
thin = 1
)
starting worker pid=3020 on localhost:11086 at 19:39:17.344
starting worker pid=3057 on localhost:11086 at 19:39:17.454
starting worker pid=3094 on localhost:11086 at 19:39:17.563
SAMPLING FOR MODEL '9afeac097877f56e1844b004ab8095d1' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 2e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 300 [ 0%] (Warmup)
Chain 1: Iteration: 30 / 300 [ 10%] (Warmup)
Chain 1: Iteration: 60 / 300 [ 20%] (Warmup)
Chain 1: Iteration: 90 / 300 [ 30%] (Warmup)
Chain 1: Iteration: 120 / 300 [ 40%] (Warmup)
Chain 1: Iteration: 150 / 300 [ 50%] (Warmup)
Chain 1: Iteration: 180 / 300 [ 60%] (Warmup)
Chain 1: Iteration: 181 / 300 [ 60%] (Sampling)
Chain 1: Iteration: 210 / 300 [ 70%] (Sampling)
Chain 1: Iteration: 240 / 300 [ 80%] (Sampling)
Chain 1: Iteration: 270 / 300 [ 90%] (Sampling)
Chain 1: Iteration: 300 / 300 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.422889 seconds (Warm-up)
Chain 1: 0.113815 seconds (Sampling)
Chain 1: 0.536704 seconds (Total)
Chain 1:
SAMPLING FOR MODEL '9afeac097877f56e1844b004ab8095d1' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 2e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 300 [ 0%] (Warmup)
Chain 2: Iteration: 30 / 300 [ 10%] (Warmup)
Chain 2: Iteration: 60 / 300 [ 20%] (Warmup)
Chain 2: Iteration: 90 / 300 [ 30%] (Warmup)
Chain 2: Iteration: 120 / 300 [ 40%] (Warmup)
Chain 2: Iteration: 150 / 300 [ 50%] (Warmup)
Chain 2: Iteration: 180 / 300 [ 60%] (Warmup)
Chain 2: Iteration: 181 / 300 [ 60%] (Sampling)
SAMPLING FOR MODEL '9afeac097877f56e1844b004ab8095d1' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 1.5e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 300 [ 0%] (Warmup)
Chain 3: Iteration: 30 / 300 [ 10%] (Warmup)
Chain 2: Iteration: 210 / 300 [ 70%] (Sampling)
Chain 2: Iteration: 240 / 300 [ 80%] (Sampling)
Chain 2: Iteration: 270 / 300 [ 90%] (Sampling)
Chain 2: Iteration: 300 / 300 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.499519 seconds (Warm-up)
Chain 2: 0.129509 seconds (Sampling)
Chain 2: 0.629028 seconds (Total)
Chain 2:
Chain 3: Iteration: 60 / 300 [ 20%] (Warmup)
Chain 3: Iteration: 90 / 300 [ 30%] (Warmup)
Chain 3: Iteration: 120 / 300 [ 40%] (Warmup)
Chain 3: Iteration: 150 / 300 [ 50%] (Warmup)
Chain 3: Iteration: 180 / 300 [ 60%] (Warmup)
Chain 3: Iteration: 181 / 300 [ 60%] (Sampling)
Chain 3: Iteration: 210 / 300 [ 70%] (Sampling)
Chain 3: Iteration: 240 / 300 [ 80%] (Sampling)
Chain 3: Iteration: 270 / 300 [ 90%] (Sampling)
Chain 3: Iteration: 300 / 300 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.498432 seconds (Warm-up)
Chain 3: 0.1909 seconds (Sampling)
Chain 3: 0.689332 seconds (Total)
Chain 3:
params_2 <- c('beta[1]', 'beta[2]', 'beta[3]', 'beta[4]', 'beta[5]', 'beta[6]', 'sigma')
traceplot(bayesion_model_2, inc_warmup = TRUE, pars = params_2)
lineal_model_2$coefficients
(Intercept) bill_length_mm bill_depth_mm flipper_length_mm sexmale
-2399.959627 -1.251644 -81.390712 38.730138 509.703109
br_coeficients(bayesion_model_2, params_2)
vars_2 <- c('bill_length_mm', 'bill_depth_mm', 'flipper_length_mm', 'sex_female','sex_male')
lm_vs_br_models_validation(
lineal_model_2,
bayesion_model_2,
params_2,
vars_2,
test_set,
test_set_cat
)
bayesion_predictor_2 <- BayesianRegressionPredictor.from(bayesion_model_2, params_2, vars_2)
plot_compare_fit(
lineal_model_2,
bayesion_predictor_2,
test_set,
test_set_cat,
label_1='Regresion Lineal',
label_2='Regresion Bayesiana'
)
A continuación vamos a agregar outliers a la variable flipper_length_mm, la cual define la longitud de la aleta del individuo medida en milÃmetros.
Para visualizar los nuevos outliers a continuación graficamos flipper_length_mm vs. body_mass_g antes y después de agregar outliers:
plot_data(train_set)
train_set_with_outliers <- train_set %>%
mutate(flipper_length_mm = ifelse(
body_mass_g > 4900 & body_mass_g < 5000,
flipper_length_mm + (flipper_length_mm * runif(1, 0.1, 0.2)),
flipper_length_mm
))
plot_data(train_set_with_outliers)
lineal_model_3 <- lm(
body_mass_g ~ bill_length_mm + bill_depth_mm + flipper_length_mm,
data = train_set_with_outliers
)
Comparemos el ajuste del modelo lineal ajustando en un dataset de train con y sin outliers:
plot_compare_fit(
lineal_model_1,
lineal_model_3,
train_set_with_outliers,
label_1='Regresión Lineal SIN outliers',
label_2='Regresión Lineal CON outliters'
)
Se aprecia que el modelo entrenado en el train set outliers esta apalancado por las observaciones atipicas.
Entrenamos una regresión lineal múltiple robusta para intentar de disminuir el efecto de los nuevos outliers.
p_load(MASS)
robust_lineal_model_3 <- rlm(
body_mass_g ~ bill_length_mm + bill_depth_mm + flipper_length_mm,
data = train_set_with_outliers
)
plot_compare_fit(
lineal_model_1,
robust_lineal_model_3,
train_set_with_outliers,
label_1='Regresión Lineal SIN outliers',
label_2='Regresión Lineal Robusta CON outliters'
)
Gráficamente no se llega a distinguir pero el modelo robusto termina ajustando mejor que el modelo lineal clásico. Esto se puede observar cuando comparamos el RMSE/MAE sobre train.
lm_vs_lm_models_validation(lineal_model_3, robust_lineal_model_3, train_set_with_outliers)
Warning in actual - predicted :
longer object length is not a multiple of shorter object length
Warning in actual - predicted :
longer object length is not a multiple of shorter object length
Warning in actual - predicted :
longer object length is not a multiple of shorter object length
Warning in actual - predicted :
longer object length is not a multiple of shorter object length
Mas alla de esto el modelo clasico sigue dando mejores resultado en test:
lm_vs_lm_models_validation(lineal_model_3, robust_lineal_model_3, test_set)
bayesion_model_3 <- stan(
model_code = "
data {
int<lower=1> obs_count;
vector<lower=1>[obs_count] x1;
vector<lower=1>[obs_count] x2;
vector<lower=1>[obs_count] x3;
vector[obs_count] y;
}
parameters {
real beta0;
real beta1;
real beta2;
real beta3;
real<lower=0> sigma;
}
model {
// Distribuciones a priori
beta0 ~ normal(0, 8000); // Intecept
beta1 ~ normal(0, 100);
beta2 ~ normal(0, 100);
beta3 ~ normal(0, 100);
sigma ~ exponential(0.1);
// Likelihood
y ~ normal(beta0 + beta1 * x1 + beta2 * x2 + beta3 * x3, sigma);
}
",
data = list(
obs_count = nrow(train_set_with_outliers),
y = colvalues(train_set_with_outliers, 'body_mass_g'),
x1 = colvalues(train_set_with_outliers, 'bill_length_mm'),
x2 = colvalues(train_set_with_outliers, 'bill_depth_mm'),
x3 = colvalues(train_set_with_outliers, 'flipper_length_mm')
),
chains = 3,
iter = 300,
warmup = 180,
thin = 1
)
starting worker pid=3296 on localhost:11086 at 19:39:48.793
starting worker pid=3333 on localhost:11086 at 19:39:48.903
starting worker pid=3370 on localhost:11086 at 19:39:49.012
SAMPLING FOR MODEL 'c0c5a3ba8590ba28e2cb449cd35b7ade' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 3.3e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.33 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 300 [ 0%] (Warmup)
Chain 1: Iteration: 30 / 300 [ 10%] (Warmup)
Chain 1: Iteration: 60 / 300 [ 20%] (Warmup)
Chain 1: Iteration: 90 / 300 [ 30%] (Warmup)
Chain 1: Iteration: 120 / 300 [ 40%] (Warmup)
SAMPLING FOR MODEL 'c0c5a3ba8590ba28e2cb449cd35b7ade' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 4.7e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.47 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 300 [ 0%] (Warmup)
Chain 2: Iteration: 30 / 300 [ 10%] (Warmup)
Chain 1: Iteration: 150 / 300 [ 50%] (Warmup)
Chain 1: Iteration: 180 / 300 [ 60%] (Warmup)
Chain 1: Iteration: 181 / 300 [ 60%] (Sampling)
Chain 2: Iteration: 60 / 300 [ 20%] (Warmup)
Chain 1: Iteration: 210 / 300 [ 70%] (Sampling)
Chain 1: Iteration: 240 / 300 [ 80%] (Sampling)
SAMPLING FOR MODEL 'c0c5a3ba8590ba28e2cb449cd35b7ade' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 4.3e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.43 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 300 [ 0%] (Warmup)
Chain 3: Iteration: 30 / 300 [ 10%] (Warmup)
Chain 1: Iteration: 270 / 300 [ 90%] (Sampling)
Chain 1: Iteration: 300 / 300 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 1.31268 seconds (Warm-up)
Chain 1: 0.434985 seconds (Sampling)
Chain 1: 1.74766 seconds (Total)
Chain 1:
Chain 2: Iteration: 90 / 300 [ 30%] (Warmup)
Chain 3: Iteration: 60 / 300 [ 20%] (Warmup)
Chain 2: Iteration: 120 / 300 [ 40%] (Warmup)
Chain 2: Iteration: 150 / 300 [ 50%] (Warmup)
Chain 2: Iteration: 180 / 300 [ 60%] (Warmup)
Chain 2: Iteration: 181 / 300 [ 60%] (Sampling)
Chain 3: Iteration: 90 / 300 [ 30%] (Warmup)
Chain 2: Iteration: 210 / 300 [ 70%] (Sampling)
Chain 2: Iteration: 240 / 300 [ 80%] (Sampling)
Chain 3: Iteration: 120 / 300 [ 40%] (Warmup)
Chain 2: Iteration: 270 / 300 [ 90%] (Sampling)
Chain 3: Iteration: 150 / 300 [ 50%] (Warmup)
Chain 2: Iteration: 300 / 300 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 1.1051 seconds (Warm-up)
Chain 2: 0.495396 seconds (Sampling)
Chain 2: 1.6005 seconds (Total)
Chain 2:
Chain 3: Iteration: 180 / 300 [ 60%] (Warmup)
Chain 3: Iteration: 181 / 300 [ 60%] (Sampling)
Chain 3: Iteration: 210 / 300 [ 70%] (Sampling)
Chain 3: Iteration: 240 / 300 [ 80%] (Sampling)
Chain 3: Iteration: 270 / 300 [ 90%] (Sampling)
Chain 3: Iteration: 300 / 300 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 1.25576 seconds (Warm-up)
Chain 3: 0.638169 seconds (Sampling)
Chain 3: 1.89393 seconds (Total)
Chain 3:
params_3 <- c('beta0', 'beta1', 'beta2', 'beta3', 'sigma')
traceplot(bayesion_model_3, pars = params_3, inc_warmup = TRUE)
lm_vs_br_coeficients(robust_lineal_model_3, bayesion_model_3, params_3)
vars_3 <- c('bill_length_mm', 'bill_depth_mm', 'flipper_length_mm')
lm_vs_br_models_validation(robust_lineal_model_3, bayesion_model_3, params_3, vars_3, test_set)
bayesion_predictor_3 <- BayesianRegressionPredictor.from(bayesion_model_3, params_3, vars_3)
plot_compare_fit(
robust_lineal_model_3,
bayesion_predictor_3,
train_set,
label_1='Regresion Lineal Robusta CON outliers',
label_2='Regresion Bayesiana CON outliers'
)
COMPLETAR
En este aso entrenamos solo con el 10% de lo datos.
train_test <- train_test_split(dataset, train_size = 0.05, shuffle = TRUE)
[1] "Train set size: 16"
[1] "Test set size: 317"
train_set_4 <- train_test[[1]]
test_set_4 <- train_test[[2]]
plot_data(train_set_4)
lineal_model_4 <- lm(
body_mass_g ~ bill_length_mm + bill_depth_mm + flipper_length_mm,
data = train_set_4
)
bayesion_model_4 <- stan(
model_code = "
data {
int<lower=1> obs_count;
vector<lower=1>[obs_count] x1;
vector<lower=1>[obs_count] x2;
vector<lower=1>[obs_count] x3;
vector[obs_count] y;
}
parameters {
real beta0;
real beta1;
real beta2;
real beta3;
real<lower=0> sigma;
}
model {
// Distribuciones a priori
beta0 ~ normal(0, 8000); // Intercept
beta1 ~ normal(0, 100);
beta2 ~ normal(0, 100);
beta3 ~ normal(0, 100);
sigma ~ exponential(0.1);
// Likelihood
y ~ normal(beta0 + beta1 * x1 + beta2 * x2 + beta3 * x3, sigma);
}
",
data = list(
obs_count = nrow(train_set_4),
y = colvalues(train_set_4, 'body_mass_g'),
x1 = colvalues(train_set_4, 'bill_length_mm'),
x2 = colvalues(train_set_4, 'bill_depth_mm'),
x3 = colvalues(train_set_4, 'flipper_length_mm')
),
chains = 3,
iter = 300,
warmup = 180,
thin = 1
)
starting worker pid=3572 on localhost:11086 at 19:40:19.482
starting worker pid=3609 on localhost:11086 at 19:40:19.591
starting worker pid=3646 on localhost:11086 at 19:40:19.700
SAMPLING FOR MODEL 'a680c5a9703a67be8bd5c21746d8ddcd' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 1.5e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 300 [ 0%] (Warmup)
Chain 1: Iteration: 30 / 300 [ 10%] (Warmup)
Chain 1: Iteration: 60 / 300 [ 20%] (Warmup)
Chain 1: Iteration: 90 / 300 [ 30%] (Warmup)
Chain 1: Iteration: 120 / 300 [ 40%] (Warmup)
Chain 1: Iteration: 150 / 300 [ 50%] (Warmup)
Chain 1: Iteration: 180 / 300 [ 60%] (Warmup)
Chain 1: Iteration: 181 / 300 [ 60%] (Sampling)
Chain 1: Iteration: 210 / 300 [ 70%] (Sampling)
Chain 1: Iteration: 240 / 300 [ 80%] (Sampling)
Chain 1: Iteration: 270 / 300 [ 90%] (Sampling)
Chain 1: Iteration: 300 / 300 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.154591 seconds (Warm-up)
Chain 1: 0.046818 seconds (Sampling)
Chain 1: 0.201409 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'a680c5a9703a67be8bd5c21746d8ddcd' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 1.2e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 300 [ 0%] (Warmup)
Chain 2: Iteration: 30 / 300 [ 10%] (Warmup)
Chain 2: Iteration: 60 / 300 [ 20%] (Warmup)
Chain 2: Iteration: 90 / 300 [ 30%] (Warmup)
Chain 2: Iteration: 120 / 300 [ 40%] (Warmup)
Chain 2: Iteration: 150 / 300 [ 50%] (Warmup)
Chain 2: Iteration: 180 / 300 [ 60%] (Warmup)
Chain 2: Iteration: 181 / 300 [ 60%] (Sampling)
Chain 2: Iteration: 210 / 300 [ 70%] (Sampling)
Chain 2: Iteration: 240 / 300 [ 80%] (Sampling)
Chain 2: Iteration: 270 / 300 [ 90%] (Sampling)
Chain 2: Iteration: 300 / 300 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.113693 seconds (Warm-up)
Chain 2: 0.063825 seconds (Sampling)
Chain 2: 0.177518 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'a680c5a9703a67be8bd5c21746d8ddcd' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 1.1e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 300 [ 0%] (Warmup)
Chain 3: Iteration: 30 / 300 [ 10%] (Warmup)
Chain 3: Iteration: 60 / 300 [ 20%] (Warmup)
Chain 3: Iteration: 90 / 300 [ 30%] (Warmup)
Chain 3: Iteration: 120 / 300 [ 40%] (Warmup)
Chain 3: Iteration: 150 / 300 [ 50%] (Warmup)
Chain 3: Iteration: 180 / 300 [ 60%] (Warmup)
Chain 3: Iteration: 181 / 300 [ 60%] (Sampling)
Chain 3: Iteration: 210 / 300 [ 70%] (Sampling)
Chain 3: Iteration: 240 / 300 [ 80%] (Sampling)
Chain 3: Iteration: 270 / 300 [ 90%] (Sampling)
Chain 3: Iteration: 300 / 300 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.126041 seconds (Warm-up)
Chain 3: 0.024386 seconds (Sampling)
Chain 3: 0.150427 seconds (Total)
Chain 3:
params_4 <- c('beta0', 'beta1', 'beta2', 'beta3', 'sigma')
traceplot(bayesion_model_4, pars = params_4, inc_warmup = TRUE)
Coeficientes de la regresión múltiple:
lineal_model_4$coefficients
(Intercept) bill_length_mm bill_depth_mm flipper_length_mm
-8874.92522 24.46866 76.75274 53.51397
Coeficientes descubiertos por la regresión múltiple bayesiana:
for(param in params_4) print(get_posterior_mean(bayesion_model_4, par=param)[4])
[1] -8223.236
[1] 23.60412
[1] 61.40233
[1] 51.73701
[1] 241.8201
vars_4 <- c('bill_length_mm', 'bill_depth_mm', 'flipper_length_mm')
lm_vs_br_models_validation(lineal_model_4, bayesion_model_4, params_4, vars_4, test_set)
bayesion_predictor_4 <- BayesianRegressionPredictor.from(bayesion_model_4, params_4, vars_4)
plot_compare_fit(
lineal_model_4,
bayesion_predictor_4,
train_set,
label_1='Regresion Lineal',
label_2='Regresion Bayesiana'
)
Definimos una distribución uniforme para el beta asociado a la variable flipper_length_mm.
starting worker pid=3848 on localhost:11086 at 19:40:49.175
starting worker pid=3885 on localhost:11086 at 19:40:49.287
starting worker pid=3922 on localhost:11086 at 19:40:49.396
SAMPLING FOR MODEL 'ec5409959af1c66bd72a3906ea8aa743' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 3.6e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.36 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 1: Iteration: 181 / 1000 [ 18%] (Sampling)
SAMPLING FOR MODEL 'ec5409959af1c66bd72a3906ea8aa743' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 4.4e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.44 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 280 / 1000 [ 28%] (Sampling)
Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
SAMPLING FOR MODEL 'ec5409959af1c66bd72a3906ea8aa743' NOW (CHAIN 3).
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: exponential_lpdf: Random variable is -0.274745, but must be >= 0! (in 'model19f65fff8b1_ec5409959af1c66bd72a3906ea8aa743' at line 21)
Chain 3:
Chain 3: Gradient evaluation took 4e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.4 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 380 / 1000 [ 38%] (Sampling)
Chain 2: Iteration: 181 / 1000 [ 18%] (Sampling)
Chain 2: Iteration: 280 / 1000 [ 28%] (Sampling)
Chain 1: Iteration: 480 / 1000 [ 48%] (Sampling)
Chain 2: Iteration: 380 / 1000 [ 38%] (Sampling)
Chain 3: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 2: Iteration: 480 / 1000 [ 48%] (Sampling)
Chain 1: Iteration: 580 / 1000 [ 58%] (Sampling)
Chain 3: Iteration: 181 / 1000 [ 18%] (Sampling)
Chain 2: Iteration: 580 / 1000 [ 58%] (Sampling)
Chain 2: Iteration: 680 / 1000 [ 68%] (Sampling)
Chain 1: Iteration: 680 / 1000 [ 68%] (Sampling)
Chain 3: Iteration: 280 / 1000 [ 28%] (Sampling)
Chain 2: Iteration: 780 / 1000 [ 78%] (Sampling)
Chain 2: Iteration: 880 / 1000 [ 88%] (Sampling)
Chain 1: Iteration: 780 / 1000 [ 78%] (Sampling)
Chain 3: Iteration: 380 / 1000 [ 38%] (Sampling)
Chain 2: Iteration: 980 / 1000 [ 98%] (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 1.01022 seconds (Warm-up)
Chain 2: 2.02517 seconds (Sampling)
Chain 2: 3.0354 seconds (Total)
Chain 2:
Chain 1: Iteration: 880 / 1000 [ 88%] (Sampling)
Chain 3: Iteration: 480 / 1000 [ 48%] (Sampling)
Chain 1: Iteration: 980 / 1000 [ 98%] (Sampling)
Chain 3: Iteration: 580 / 1000 [ 58%] (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 1.10362 seconds (Warm-up)
Chain 1: 3.90271 seconds (Sampling)
Chain 1: 5.00633 seconds (Total)
Chain 1:
Chain 3: Iteration: 680 / 1000 [ 68%] (Sampling)
Chain 3: Iteration: 780 / 1000 [ 78%] (Sampling)
Chain 3: Iteration: 880 / 1000 [ 88%] (Sampling)
Chain 3: Iteration: 980 / 1000 [ 98%] (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 1.14468 seconds (Warm-up)
Chain 3: 3.67208 seconds (Sampling)
Chain 3: 4.81675 seconds (Total)
Chain 3:
br_vs_br_coeficients(bayesion_model_1, bayesion_model_5, params_5)
lm_vs_br_models_validation(lineal_model_1, bayesion_model_1, params_1, vars_1, test_set)
vars_5 <- c('bill_length_mm', 'bill_depth_mm', 'flipper_length_mm')
lm_vs_br_models_validation(lineal_model_1, bayesion_model_5, params_5, vars_5, test_set)
bayesion_predictor_5 <- BayesianRegressionPredictor.from(bayesion_model_5, params_5, vars_5)
plot_compare_fit(
bayesion_predictor_1,
bayesion_predictor_5,
train_set,
label_1='Regresion Bayesiana con dist informativa',
label_2='Regresion Bayesiana con dist menos informativa'
)
COMPLETAR